%%算法参数
    %%响度loudness，脉冲pulse，速度velosity,位置position
%音强衰减系数alpha，脉冲增加系数gama，频率最大值fmax，频率最小值fmin
%%初始化种群，计算适应度，找最优个体
%计算平均响度，   
%迭代更新
%迭代次数大于threshold。被其他干扰，最优值干扰或随机干扰，产生随机数看是否在最优值附近产生干扰
%如干扰后适应度较大则粒子更新否则不更新

clc;clear
%% 下载数据
% 加工数据包括加工时间，   加工机器，         机器数，    各机器权重，   工件数，   各工件对应的工序数

global  operation_time operation_machine num_machine num_job num_op;
load Mk01 operation_time operation_machine num_machine  num_job num_op 

%% 基本参数
sizepop=50;       % 种群数量，一般20-40
MAXGEN=200;   % 最大迭代次数
loudness=0.8;            % 初始响度，可更大，搜索范围
pulse=0.4;           % 初始脉冲，可更大，搜索频度
alpha=0.99;     %音强衰减系数
gamma=0.1;      % 脉冲增加系数
Frequency_min=0;     % 频率最小
Frequency_max=2;     % 频率最大
gen=0;            % I初始迭代次数
n=0.05;    %混沌映射概率
total_op_num=sum(num_op);          %搜索维度---工序总数
%% ===========================种群初始化=======================================
trace = zeros(2,MAXGEN);
frequency=zeros(sizepop,1);   %频率20*1
velocities=rand(sizepop,total_op_num);     
position=rand(sizepop,total_op_num)*2-1; %位置-1~1之间
FF=initialization();
[Z_new,~,~,~,~]=fitness(chroms_new,FF);
positions_new=zeros(sizepop,total_op_num)   ;%更新时新位置

% 找最优
[fmin,I]=min(Fitness);%选最小
best=position(I,:);

%% =================================迭代过程==============================

while (gen<MAXGEN)
   % 计算响度和脉冲发散率
   pulse_r=pulse*(1.0-exp(-gamma*gen));
   loudness=alpha*loudness;  %每个个体响度相同
  % 更新个体
   for i=1:sizepop   %对于个体i
       frequency(i)=Frequency_min+(Frequency_max-Frequency_min)*rand;%频率随机产生，受最优值影响速度程度0~2间波动
       velocities(i,:)=velocities(i,:)+(position(i,:)-best)*frequency(i);   %朝着最优值走去
       positions_new(i,:)=position(i,:)+velocities(i,:);
  
       % 迭代次数越大，被最优值影响概率越大
   if rand<pulse_r
       positions_new(i,:)=best+0.1*randn(1,total_op_num)*loudness;
   end  %产生均值为0，方差为1，标准差为1的正态随机数

   end %%个体产生新解结束
   
   % 检查新位置是否在范围内
   [positions_new,velocities]=bounds(positions_new,velocities,sizepop,total_op_num);
   % 计算新位置适应度
   [Z_new,~,~,~,~]=fitness(chroms_new,FF);
   % 新位置较好并且不太远则更新位置
   for i=1:sizepop
    if ((Fitness_new(i)<=Fitness(i)) && (rand<loudness))
       position(i,:)=positions_new(i,:);
       Fitness(i)=Fitness_new(i);
    end
   % 更新最优
    if Fitness_new(i)<=fmin
       best=positions_new(i,:);
       fmin=Fitness_new(i);
    end
   end
  gen=gen+1 ;  % 迭代次数+1
   fprintf('当前迭代次数：'),disp(gen);
  % 记录每代的最优适应度与平均适应度
    trace(1, gen)=fmin;       
    trace(2, gen)=mean(Fitness);  
    % 更新全局最优适应度
    if gen==1 || MinVal>trace(1,gen)
        MinVal=trace(1,gen);
    end
end  % 迭代结束
%% ===============================输出结果==================================
%输出最好位置
disp(['Best =',num2str(best),'        fmin=',num2str(fmin)]);

%% ============================输出结果=============================
%% 输出最优适应度
fprintf('最优适应度：'),disp(MinVal)
%% 描绘解的变化
figure(1)
plot(trace(1,:));
hold on;
plot(trace(2,:),'-.');grid;
legend('解的变化','种群均值的变化');
%% 显示最优解
[~,~,Pvals,chroms]=fitness(chroms_new,FF);
Pval1=Pvals{1,1};
figure(2);
for i=1:total_op_num
    mText=chroms(1,total_op_num+i);
    b=chroms(1,i);
    x1=Pval1(1,i);
    x2=Pval1(2,i); 
    y1=mText-0.2;
    y2=mText;
    hold on; 
    fill([x1,x2,x2,x1],[y1,y1,y2,y2],[1-1/b,1/b,b/num_job]);
    text((x1+x2)/2,mText-0.1,num2str(b));
end